Q3:

image-2.png

image.png

In [2]:
pip install mermaid-py
Requirement already satisfied: mermaid-py in /home/nicai_zht/miniconda3/envs/sc-python/lib/python3.12/site-packages (0.8.0)
Requirement already satisfied: requests<3.0.0,>=2.31.0 in /home/nicai_zht/miniconda3/envs/sc-python/lib/python3.12/site-packages (from mermaid-py) (2.32.3)
Requirement already satisfied: charset-normalizer<4,>=2 in /home/nicai_zht/miniconda3/envs/sc-python/lib/python3.12/site-packages (from requests<3.0.0,>=2.31.0->mermaid-py) (3.4.2)
Requirement already satisfied: idna<4,>=2.5 in /home/nicai_zht/miniconda3/envs/sc-python/lib/python3.12/site-packages (from requests<3.0.0,>=2.31.0->mermaid-py) (3.10)
Requirement already satisfied: urllib3<3,>=1.21.1 in /home/nicai_zht/miniconda3/envs/sc-python/lib/python3.12/site-packages (from requests<3.0.0,>=2.31.0->mermaid-py) (2.4.0)
Requirement already satisfied: certifi>=2017.4.17 in /home/nicai_zht/miniconda3/envs/sc-python/lib/python3.12/site-packages (from requests<3.0.0,>=2.31.0->mermaid-py) (2025.4.26)
Note: you may need to restart the kernel to use updated packages.
In [3]:
import mermaid as md

Z博士正在计划一个项目,研究宿主人类基因组和肠道微生物群对治疗性食物干预的影响。

一个20万患者的队列将被收集起来。

在项目的第一阶段,200名患者将在3个月的治疗性食物干预之前和之后进行基因组测序(覆盖30倍),以及RNA-seq数据(每个约30Gbps)。

在3个月的治疗性食物干预之前和之后,所有这100名患者还将进行肠道微生物组测序(每个样本至少10GB)。

包括这200例患者,总共2,000例患者将在3个月的治疗性食物干预之前和之后进行肠道微生物的16S rRNA基因测序(每个样本至少30Mbps)。

请给出一个涵盖测序和数据分析阶段的策略。你需要展示这个策略(用图表),并给出你对这个策略的时间和成本的估计。例如,你可以选择不同的测序平台或任何你知道的测序平台。

对于16S rRNA测序,500个样本可以合并到一次illumine测序中。

对于成本估算,请给出合理的估算。例如,您可以选择云计算或上海交通大学的Pi2.0超级计算机来进行数据分析。您可以从任何计算资源提供商,如cn.aliyun.com或www.huaweicloud.com找到当前的市场价格。如果你使用Pi2.0,价格可以设定为~0.02¥/(核心小时)计算和~50¥/(TB年)存储。

如果Z博士想在上述项目的第二阶段做100倍规模的测序,在时间和金钱方面的成本是多少?根据你最好的知识来回答这些问题。

1,测序阶段策略

主要是各个测序平台的比较,本质上就是benchmark

(1)基因组测序 :WGS、WES

DNA测序平台比较:

image.png

image-3.png

ref文献:Performance assessment of DNA sequencing platforms in the ABRF Next-Generation Sequencing Study

image-2.png

image.png

image-3.png

image-4.png

关于该测序平台的官网指标:

参考:https://en.mgi-tech.com/products/instruments_info/5/

image.png

image-2.png

image-3.png

测序平台主要是国内华大BGI的以及国外的Illumina:

对于 200 例患者基因组测序(30x 覆盖度),可选择华大的 DNBSEQ-T7 平台,其具有高通量和相对较低的成本,能在较短时间内完成大量样本的测序。

根据样本数量和覆盖度要求,估算所需的测序芯片数量和时间。

200 例患者每例约需 90Gb 数据(基于人类基因组大小约 3Gb,30x 覆盖度),总共约需 18Tb 数据(粗略估计,1024转换略)

根据官网内容,DNBSEQ-T7测序仪的性能参数如下:

  • 单次运行最大数据量:7 Tb(使用PE150测序模式,基于V3.0试剂套装)。
  • 运行时间:22小时至24小时(PE150模式)。

时间估计¶

  1. 总数据需求:200例患者 × 90 Gb/例 = 18 Tb。
  2. 单次运行数据量:7 Tb。
  3. 所需运行次数:18 Tb ÷ 7 Tb/次 ≈ 2.57次,向上取整为3次。
  4. 总运行时间:3次 × 24小时/次 = 72小时(假设每次运行时间均为24小时,且不考虑中间停机时间)。——》当然我们可以并行,此处计算的是“总运行机时”

(2)RNA-seq 转录组测序

转录组测序 :对于 200 例患者的 RNA-seq 数据(每例约 30Gbps),同样可采用华大的 DNBSEQ-T7 平台(也可以做转录组测序)。

时间估计¶

  1. 总数据需求:200例患者 × 30 Gb/例 = 6 Tb。
  2. 单次运行数据量:7 Tb。
  3. 所需运行次数:6 Tb ÷ 7 Tb/次 ,向上取整为1次。
  4. 总运行时间:1次 × 24小时/次 = 24小时(假设每次运行时间均为24小时,且不考虑中间停机时间)。

(3)肠道微生物组测序:宏基因组

image.png

image-2.png

对于宏基因组测序的数据,也同样可以使用该平台,保持一致(也算是个批次效应),

时间估计¶

  1. 总数据需求:100例患者 × 10 Gb/例 = 1 Tb。(题目中说的是100例患者,我觉得可能说错了,应该还还是200名患者,毕竟干预前后应该都要做一次测序收集数据,所以人数应该都是统一的,其他的就是统一的)
  2. 单次运行数据量:7 Tb。
  3. 所需运行次数:1 Tb ÷ 7 Tb/次 ,向上取整为1次。(如果是200名患者,实际上还是可以在一次运行中完成,所以时间成本还是1天)
  4. 总运行时间:1次 × 24小时/次 = 24小时(假设每次运行时间均为24小时,且不考虑中间停机时间)。

(4) 16s rRNA测序:

如果还是使用华大的平台的话,可以使用CycloneSEQ-WT02,这个是纳米孔测序,属于是三代测序范畴了,前面几个组学使用的都是成熟的二代测序平台;

image.png

参考:https://en.mgitech.cn/Home/Products/instruments_info/id/65.html

image-2.png

image-3.png

如果还是使用二代测序平台的话,华大的暂时没找到,可以使用Illumina测序平,

参考:https://www.illumina.com.cn/areas-of-interest/microbiology/microbial-sequencing-methods/16s-rrna-sequencing.html

image.png

可以使用NextSeq 2000测序平台

对于16s rRNA扩增子测序的数据,我们可以使用华大的三代测序平台:

时间估计¶

  1. 总数据需求:2000例患者 × 0.03 Gb/例 = 0.06 Tb = 60 Gb。
  2. 单次运行数据量:50 Gb。我们主要看最大数据吞吐量的运行时间效率,比如说官网这里说每个flow cell可以输出10Mb/min,主要是看100Mb/10mins的水平
  3. 所需运行次数:60 Gb ÷ 10 Mb/min = 100 h
  4. 总运行时间:该机器最大运行时间为72h,也就是3天,所以我们需要分开数据运行,串行2次或者并行,总之100h,大约5天左右,我们计算的是总运行机时

综上:¶

总时间成本:测序要在3个月前后都测一遍,所以是2次成本;

单次成本的话:

全基因组测序72h,转录组测序24h,肠道微生物宏基因组24h(无论题目中说是这100例还是200例),扩增子测序100h,单次总共是220h,也就是10天机时

总共前后3个月的话,就是440h,也就是20天机时

如果是开销成本的话:

没找到明确的测序开销成本,但是测序仪器的成本倒是能够查到,无法估计,用测序平台的仪器,但是送测任务会依据承包商上下波动

对于二代二序平台:也就是DNBSEQ-T7

image-2.png

成本上是每个基因组150美元,

按照现在的美元-人民币汇率:取1美元=7.19人民币

image-3.png

我们按照数据量来看,如果每个基因组3G就是150美元(50美元/G)的话,那么前面使用二代测序平台的数据量一共有:(200x30x3+200x30+100x10+2000x0.03)x 2 = 50120 G, 50T左右的数据量

50120*50 = 2506000美元

参考:http://completegenomics.com/products/sequencing-platforms/dnbseq-t7/

对于三代测平台,也就是

参考:https://mgi-whatif.com/eoypromo2024

image.png

我们可以使用前面next-seq 2000的价格来估计,参考https://genohub.com/ngs-sequencer/7/illumina-nextseq-1000/2000/

我们不买仪器,就送测,粗略估计成本就在测序成本上下:

image-4.png

大概79.69美元/G数据,扩增子测序数据量在60Gx2上,所以是9562.8美元

综上粗略估计2515562.8美元,约合18086896.531999998人民币,即1800万人民币左右

2,数据分析阶段

测序阶段的成本估计:

时间成本还是消费成本,实际上就是看测序平台(测序仪器的数据吞吐量决定运行时间,运行时间决定消费成本);

数据分析阶段的成本估计:

这里应该是在大数据规模下,进行1个时间复杂度的分析,涉及到超算计算任务运行收费;另外还有空间复杂度的分析,因为冷存储也要收费,冷存储的时间成本看实验周期也就是项目周期跨度;

(1)基因组测序:WGS/WES

In [6]:
%%mermaidjs # with --img flag in case your notebook doesn't render HTML
flowchart LR
    A-->B
    B-->C
No description has been provided for this image

1.png

In [8]:
%%mermaidjs
graph TD
    A[原始FASTQ文件] --> B[Read比对]
    B --> C[排序]
    C --> D[去重复]
    D --> E[Indel区域重比对]
    E --> F[碱基质量值重校正<br/><i>BQSR</i>]
    F --> G[变异检测]
    G --> H[变异质控和过滤<br/><i>VQSR</i>]
    H --> I[QC<br/>质量评估]

原始FASTQ文件

Read比对

排序

去重复

Indel区域重比对

碱基质量值重校正
BQSR

变异检测

变异质控和过滤
VQSR

QC
质量评估

WGS中变异检测,其实用到最多的工具大概就是bcftools和GATK套件了

bcftools:https://samtools.github.io/bcftools/

GATK:https://gatk.broadinstitute.org/hc/en-us

原始FASTQ:测序仪器下机的原始数据文件

Read比对:将测序片段比对到参考基因组(常用BWA/Bowtie2)

排序:按染色体坐标排序(常用samtools sort)

去重复:移除PCR重复序列(常用Picard MarkDuplicates/samtools也有rmdup类似子命令)

Indel区域重比对:优化插入/缺失区域比对(GATK RealignerTargetCreator)

BQSR:碱基质量值重校正,即Base (Quality Score) Recalibration (GATK BaseRecalibrator)

变异检测:识别SNP/InDel(GATK HaplotypeCaller/bcftools mpileup+call)

VQSR:基于机器学习模型的变异过滤,即Variant Quality Score Recalibration,VQSR 不会对原来的质量分数值进行校正,而是将 QUAL 列中的分值没有考虑的各项与变异相关的属性纳入其中,重新计算一个新的质量分数 VQSLOD (GATK VariantRecalibrator)

QC:最终结果的质量评估(Qualimap/MultiQC/或者是常用的fastp那些)

(2)bulk RNA-seq

In [ ]:
%%mermaidjs
graph TD
    A[Quality Control] -->|FastQC, fastp| B[Reads Mapping]
    B -->|Unspliced:<br>BWA, Bowtie, Bowtie2| C[SAM/BAM Files]
    B -->|Spliced:<br>STAR, TopHat, Hisat2| C
    C -->|Viewed by IGV| D[Expression Quantification]
    D -->|Summarize counts:<br>featureCounts| E[DE Testing]
    D -->|FPKM/RPKM:<br>HTseq| E
    E -->|DESeq2, edgeR| F[Functional Interpretation]
    E -->|Cuffdifff| F
    F -->|Function enrichment| G[Biological Insights]
    F -->|Infer networks| G
    F -->|Data integration| G

FastQC, fastp

Unspliced:
BWA, Bowtie, Bowtie2

Spliced:
STAR, TopHat, Hisat2

Viewed by IGV

Summarize counts:
featureCounts

FPKM/RPKM:
HTseq

DESeq2, edgeR

Cuffdifff

Function enrichment

Infer networks

Data integration

Quality Control

Reads Mapping

SAM/BAM Files

Expression Quantification

DE Testing

Functional Interpretation

Biological Insights

In [4]:
%%mermaidjs
graph TD
    A[Quality Control<br/>FastQC, fastp] -->|clean reads| B(Reads Mapping)
    B -->|Unspliced alignment<br/>BWA, Bowtie, Bowtie2| C{SAM/BAM Files}
    B -->|Spliced alignment<br/>STAR, TopHat, Hisat2| C
    C -->|IGV visualization| D[[Expression Quantification]]
    D -->|Count summarization<br/>featureCounts| E[(DE Testing)]
    D -->|FPKM/RPKM calculation<br/>HTseq| E
    E -->|Statistical analysis<br/>DESeq2, edgeR| F>Functional Interpretation ]
    E -->|Alternative method<br/>Cuffdiff| F
    F -->|Enrichment analysis| G((Biological Insights))
    F -->|Network inference| G
    F -->|Multi-omics integration| G
    
    classDef qcClass fill:#e1f5fe,stroke:#01579b,stroke-width:3px
    classDef mappingClass fill:#f3e5f5,stroke:#4a148c,stroke-width:3px
    classDef quantClass fill:#e8f5e8,stroke:#1b5e20,stroke-width:3px
    classDef deClass fill:#fff3e0,stroke:#e65100,stroke-width:3px
    classDef funcClass fill:#fce4ec,stroke:#880e4f,stroke-width:3px
    classDef insightClass fill:#f1f8e9,stroke:#33691e,stroke-width:3px
    
    class A qcClass
    class B mappingClass
    class C,D quantClass
    class E deClass
    class F funcClass
    class G insightClass

clean reads

Unspliced alignment
BWA, Bowtie, Bowtie2

Spliced alignment
STAR, TopHat, Hisat2

IGV visualization

Count summarization
featureCounts

FPKM/RPKM calculation
HTseq

Statistical analysis
DESeq2, edgeR

Alternative method
Cuffdiff

Enrichment analysis

Network inference

Multi-omics integration

Quality Control
FastQC, fastp

Reads Mapping

SAM/BAM Files

Expression Quantification

DE Testing

Functional Interpretation

Biological Insights

(3)肠道微生物组测序/宏基因组测序

image.png

参考:http://abc.cbi.pku.edu.cn/seminar/PKU22_G07.pdf

image.png

参考:https://biosys.bgi.com/help/metagenomics/analytic/1-1-general-intro.html#%E7%94%9F%E4%BF%A1%E5%88%86%E6%9E%90%E6%B5%81%E7%A8%8B

  • read 质控和去宿主:使用 SOAPnuke 软件对原始数据进行过滤质控,使用 Bowtie2 比对宿主序列并去除比中的序列,生成 Clean Data。
  • 功能注释:先使用 MEGAHIT 对完成质控和去宿主的序列进行基于 k-mer 的组装,生成 contigs,再用 MetaGeneMark 软件,预测 contigs 中的基因序列。使用 CD-HIT 软件,对得到的基因进行去冗余,使用 Salmon 软件统计各个基因相对丰度表,再使用 DIAMOND 或 RGI 将非冗余基因比对到 eggNOG、KEGG、BacMet、CARD、COG、CAZy、Swiss-prot 等数据库,完成基因功能注释。
  • 物种注释:使用 Kraken2 和自建数据库(筛选 NCBI NT 数据库或 UHGG 数据库)比对来计算样品中所含有物种的序列数,再用 Bracken2 来对样品中物种的实际丰度进行估计,完成物种注释。
  • 基于基因丰度表、物种丰度表和功能丰度表,可以对基因、物种和功能分布情况进行可视化,计算其 Alpha 多样性、Beta 多样性,进行物种或功能丰度聚类分析,PCA、PLSDA、PCoA 和 NMDS 降维分析,样品聚类分析,使用 Wilcoxon/Kruscal,T 检验/方差分析、LEfSe 分析以及 KEGG 通路富集分析(Reportscore 方法)挖掘样品之间的物种组成和功能组成差异。
  • 结合环境因子、病理指标或特殊表型进行深入关联研究,能够为进一步深入研究和利用样品的物种和功能提供理论依据。

(4)16s rRNA测序

image.png

image.png

参考:https://www.omicstudio.cn/ueditor/php/upload/file/20200430/59416S%20rDNA%E6%B5%8B%E5%BA%8F%E5%88%86%E6%9E%90%E6%8A%A5%E5%91%8A%E4%B8%AD%E6%96%87%E5%B7%A5%E4%B8%9A%E7%89%88.pdf

计算核时消费需要实际数据操作,做小规模数据集测试,再用时间复杂度评估时间,另外涉及到CPU 资源按核时收费,GPU 资源按卡时收费;

交我算平台集群总费用为CPU,GPU和存储费用之和;

存储费用看数据规模,见前计算出来粗略50120G,即粗略50.12T(其中1024各种转换略),每TB每年50元;

总体来说3个月前后测序,前面说机时20天运行机时,粗略估计4个月stage1,当然数据还要分析、写文投递改稿,时间会更长;

4个月的话,50*50.12/3=835.3333333333334,

stage2的话,再扩大100倍,机时如果并行的话,只是存储83533.33333333334,也得8.35万左右

实际时间会是季度的n倍,按年来算

In [7]:
50*50.12/3*100
Out[7]:
83533.33333333334